#include <iostream>
#include <stdexcept>
#include <vector>
#include <cstdlib>
#include <ginac/ginac.h>
#include "sector_decomposition/sector_decomposition.h"

int main(int argc, char **argv)
{
  // This is an advanced example.
  // The program calculates individual sectors of the triple box.
  // Usually the available memory is not sufficient to calculate the 
  // triple box integral in one run.
  //
  // The program is used as follows:
  //  triplebox n1 n2 n3 n4 n5 ...
  // where n1, n2, n3, ... are numbers.
  // This will select from the first set of the decomposition only the
  // n1-th sector, from the subsequent set of decomposition only the
  // n2-th sector is selected and so on.
  //
  // An example would be
  //  triplebox 0 0 0 0 0 0 0
  // This selects from all sets always the first sub-sector. 
  // The tree of subsectors has at this point a depth of 7,
  // therefore in this example exactly one subsector integral without
  // further decomposition is calculated.

  try{

    using namespace sector_decomposition;
    using namespace GiNaC;

    // --------------------------------------------------------------
    int number_param = argc - 1 ;
    std::vector<int> sector_selection(number_param);

    for (int j=0; j<number_param; j++) sector_selection[j] = atoi(argv[j+1]);

    std::cout << "triplebox";
    for (int j=0; j<number_param; j++) std::cout << " " << sector_selection[j];
    std::cout << std::endl;
    std::cout << std::endl;

    // --------------------------------------------------------------
    int verbose_level = 1;

    CHOICE_STRATEGY = STRATEGY_X;

    // --------------------------------------------------------------
    monte_carlo_parameters mc_parameters = monte_carlo_parameters( 5, 15, 10000, 100000 );

    // --------------------------------------------------------------

    symbol eps("eps");
    int order = 0;

    int n     = 10;
    int loops = 3;

    symbol minus_sinv("minus_sinv"), minus_tinv("minus_tinv");

    std::vector<ex> parameters;
    symbol x1("x1"), x2("x2"), x3("x3"), x4("x4"), x5("x5"), x6("x6"), x7("x7"), x8("x8"), x9("x9"), x0("x0");
    parameters.push_back(x1);
    parameters.push_back(x2);
    parameters.push_back(x3);
    parameters.push_back(x4);
    parameters.push_back(x5);
    parameters.push_back(x6);
    parameters.push_back(x7);
    parameters.push_back(x8);
    parameters.push_back(x9);
    parameters.push_back(x0);

    integration_data global_data = integration_data(parameters, eps, order);

    std::vector<exponent> nu_minus_1(n,exponent(0,0));

    ex U = (x5+x7)*(x1+x2+x3+x4)*(x6+x8+x9+x0)
      + (x3+x6)*(x1+x2+x4)*(x8+x9+x0)
      + x3*x6*(x1+x2+x4+x8+x9+x0);

    ex F = 
      (
       x2*x4*( (x3+x5+x7)*(x8+x9+x0) + x6*(x3+x5+x7+x8+x9+x0) )
       + x5*x7*(x1+x2+x3+x4)*(x6+x8+x9+x0)
       + x8*x0*( (x1+x2+x4)*(x5+x6+x7) + x3*(x1+x2+x4+x5+x6+x7) )
       + (x2*x7+x4*x5)*x3*(x6+x8+x9+x0)
       + (x5*x0+x7*x8)*x6*(x1+x2+x3+x4)
       + (x2*x0+x4*x8)*x3*x6
       )*minus_sinv
      +x1*x3*x6*x9*minus_tinv;

    F = F.subs( minus_sinv == 1 );
    F = F.subs( minus_tinv == 1 );

    std::vector<ex> poly_list;
    poly_list.push_back(U);
    poly_list.push_back(F);

    std::vector<exponent> c(poly_list.size());
    c[0] = exponent( n-(loops+1)*2, loops+1 );
    c[1] = exponent( -n+2*loops, -loops );
    for (int k=0; k<n; k++) 
      {
	c[0].sum_up(nu_minus_1[k]);
	c[1].subtract_off(nu_minus_1[k]);
      }

    integrand my_integrand = integrand(nu_minus_1, poly_list, c);

    // --------------------------------------------------------------

    basic_sector_decomposition T;

    time_t time_start;
    time_t time_end;

    if ( verbose_level>1 )
      {
	std::cout << "Input integrand:" << std::endl;
	std::cout << my_integrand.get_integrand(global_data) << std::endl;
	std::cout << std::endl;
      }

    // --------------------------------------------------------------
    // test for homogeneous polynomials
    time(&time_start);
    if ( verbose_level>0 ) std::cout << "Homogeneous polynomials:" << std::flush;

    integrand integrand_homogeneous;

    T.generate_homogenous_polynomials(global_data, my_integrand, integrand_homogeneous);

    time(&time_end);
    if ( verbose_level>0 ) std::cout << " Done. Time in seconds: " << difftime(time_end,time_start) << std::endl;

    if ( verbose_level>1 )
      {
	std::cout << integrand_homogeneous.get_integrand(global_data) << std::endl;
	std::cout << std::endl;
      }

    // --------------------------------------------------------------
    // primary sectors
    time(&time_start);
    if ( verbose_level>0 ) std::cout << "Primary sectors:        " << std::flush;

    std::vector<integrand> primary_sectors;

    T.generate_primary_sectors(global_data, integrand_homogeneous, primary_sectors);

    time(&time_end);
    if ( verbose_level>0 ) std::cout << " Done. Time in seconds: " << difftime(time_end,time_start) << std::endl;
    if ( verbose_level>0 ) std::cout << " " << primary_sectors.size() << " primary sectors." << std::endl;

    if ( verbose_level>1 )
      {
	for (int k=0; k<primary_sectors.size(); k++)
	  {
	    std::cout << k << ": " << (primary_sectors[k]).get_integrand(global_data) << std::endl;
	  }
	std::cout << std::endl;
      }

    // --------------------------------------------------------------

    integrand local_integrand;
    if ( number_param > 0 )
      {
	local_integrand = primary_sectors[ sector_selection[0] ];
	primary_sectors.clear();
	primary_sectors.push_back(local_integrand);
      }

    // --------------------------------------------------------------

    for ( int depth=1; depth<number_param; depth++)
      {
	local_integrand = primary_sectors.back();
	primary_sectors.pop_back();

	bool flag_continue = true;

	size_t n = local_integrand.index.size();
	size_t r = local_integrand.polynomial_list.size();

	// ----------------------------------------------------
	// one variable factorizes
	do 
	  {
	    flag_continue = true;
	    multi_index_ordered index_k(n,1);
	    for( index_k.init(); !index_k.overflow(); index_k++)
	      {
		for (int j=0; j<r; j++) 
		  {
		    if ( local_integrand.check_zero(j, index_k, global_data) ) 
		      {
			// can factor the k'th variable from the j'the polynomial
			size_t k = index_k[0];
			local_integrand.polynomial_list[j] = normal(local_integrand.polynomial_list[j]
								    /global_data.list_feynman_parameter[local_integrand.index[k]]);

			(local_integrand.weight[k]).sum_up( local_integrand.poly_exponent[j] );

			flag_continue = false;
		      } // check_zero
		  } // j=0; j<r;
	      } // loop over index_k
	  } while ( !flag_continue );

	  flag_continue = true;
	  // ----------------------------------------------------
	for (int j=0; j<r; j++) 
	  {
	    for (int m=2; m<=n; m++)
	      {
		multi_index_ordered index_set(n,m);
		for( index_set.init(); !index_set.overflow(); index_set++)
		  {
		    if ( flag_continue )
		      {
			if ( local_integrand.check_zero(j, index_set, global_data) )
			  {
			    std::cout << "S: " ;
			    for (int k=0; k<index_set.size(); k++) std::cout << global_data.list_feynman_parameter[local_integrand.index[index_set[k]]] << " ";
			    std::cout << std::endl;

			    // create m new subsectors
			    for (int k=0; k<m; k++)
			      {
				integrand sub_integrand = local_integrand;

				sub_integrand.create_subsector(k, index_set, global_data);

				primary_sectors.push_back(sub_integrand);
			      } // k=0; k<m;

			    flag_continue = false;
			  } // check_zero
		      } // flag_continue
		  } // loop over index_set
	      } // m=2; m<=n;
	  } // j=0; j<r;

	local_integrand = primary_sectors[ sector_selection[depth] ];
	primary_sectors.clear();
	primary_sectors.push_back(local_integrand);
      }

    // --------------------------------------------------------------
    // iteration
    time(&time_start);
    if ( verbose_level>0 ) std::cout << "Iterated sectors:       " << std::flush;

    std::vector<integrand> sectors_out;

    T.iterate_sector_decomposition(global_data, primary_sectors, sectors_out);

    time(&time_end);
    if ( verbose_level>0 ) std::cout << " Done. Time in seconds: " << difftime(time_end,time_start) << std::endl;
    if ( verbose_level>0 ) std::cout << " " << sectors_out.size() << " sub-sectors." << std::endl;

    if ( verbose_level>1 )
      {
	for (int k=0; k<sectors_out.size(); k++)
	  {
	    std::cout << k << ": " << (sectors_out[k]).get_integrand(global_data) << std::endl;
	  }
	std::cout << std::endl;
      }

    // --------------------------------------------------------------
    // subtraction
    time(&time_start);
    if ( verbose_level>0 ) std::cout << "Subtraction terms:      " << std::flush;

    std::vector<GiNaC::ex> subtracted_terms;

    T.generate_subtraction_terms(global_data, sectors_out, subtracted_terms);

    time(&time_end);
    if ( verbose_level>0 ) std::cout << " Done. Time in seconds: " << difftime(time_end,time_start) << std::endl;
    if ( verbose_level>0 ) std::cout << " " << subtracted_terms.size() << " terms." << std::endl;

    if ( verbose_level>1 )
      {
	for (int k=0; k<subtracted_terms.size(); k++)
	  {
	    std::cout << k << ": " << subtracted_terms[k] << std::endl;
	  }
	std::cout << std::endl;
      }

    // --------------------------------------------------------------
    // expansion
    time(&time_start);
    if ( verbose_level>0 ) std::cout << "Expansion:              " << std::flush;

    GiNaC::ex expansion_out;

    T.expand_in_epsilon(global_data, subtracted_terms, expansion_out);

    time(&time_end);
    if ( verbose_level>0 ) std::cout << " Done. Time in seconds: " << difftime(time_end,time_start) << std::endl;

    if ( verbose_level>1 )
      {
	std::cout << expansion_out << std::endl;
	std::cout << std::endl;
      }

    // --------------------------------------------------------------
    // numerical evaluation
    time(&time_start);
    if ( verbose_level>0 ) std::cout << "Numerical integration:  " << std::flush;

    monte_carlo_result res;

    T.numerical_integration(global_data, mc_parameters, expansion_out, res, verbose_level);

    time(&time_end);
    if ( verbose_level>0 ) std::cout << " Done. Time in seconds: " << difftime(time_end,time_start) << std::endl;

    if ( verbose_level>0 )
      {
	std::cout << "Order " << pow(global_data.epsilon,global_data.order) << ": " 
		  << res.get_mean() << " +/- " << res.get_error() << std::endl;
      }

    // --------------------------------------------------------------

  } catch (std::exception &e)
    {
      std::cout << "Exception : " << e.what() << std::endl;
    }

  return 0;
}

